%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Decomposition.m   (Script to generate Table 10)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%--------------------------------------------------------------------------
% 1) Declare parameters
%--------------------------------------------------------------------------


global alphah alphal mh  ml Q  nu
global zeta lambda eta delta beta delta_p
global zbar chi_a chi_b rho_a rho_b gamma_a gamma_b mu iota


%--------------------------------------------------------------------------
% 2) Compute the baseline moments (M0) for 1981–85
%--------------------------------------------------------------------------

M0 = struct();
M0.basic_share=basic_share_initial;
M0.scope_avg = omega_avg_initial;     % omega_avg for 1981–85
M0.RD_H    = rd_Y_h_initial;           % rd_Y_h for 1981–85
M0.RD_L    = rd_Y_l_initial;           % rd_Y_l for 1981–85
M0.trans   = trans_share_initial;      % trans_share for 1981–85
M0.growth  = g_avg_initial^(zeta/(zeta+lambda)) - 1;  % annual growth for 1981–85


%%% —————————————————————————————————————————————
%%% 3) Loop over key parameters and generate moments
%%% —————————————————————————————————————————————

deviations = zeros(6, 6);

for i = 1 : 6
    if i==1
        phi=phi_ending; theta=theta_ending; ss=ss_initial; mu=mu_initial; iota=iota_initial; gamma_a=gamma_a_initial; gamma_b=gamma_b_initial; rho_a=rho_a_initial; rho_b=rho_b_initial; chi_b=chi_b_initial;
    elseif i==2
        phi=phi_ending; theta=theta_initial; ss=ss_initial; mu=mu_initial; iota=iota_initial; gamma_a=gamma_a_initial; gamma_b=gamma_b_initial; rho_a=rho_a_initial; rho_b=rho_b_initial; chi_b=chi_b_initial;
    elseif i==3
        phi=phi_initial; theta=theta_ending; ss=ss_initial; mu=mu_initial; iota=iota_initial; gamma_a=gamma_a_initial; gamma_b=gamma_b_initial; rho_a=rho_a_initial; rho_b=rho_b_initial; chi_b=chi_b_initial;
    elseif i==4
        phi=phi_initial; theta=theta_initial; ss=ss_ending; mu=mu_initial; iota=iota_initial; gamma_a=gamma_a_initial; gamma_b=gamma_b_initial; rho_a=rho_a_initial; rho_b=rho_b_initial; chi_b=chi_b_initial;
    elseif i==5
        phi=phi_initial; theta=theta_initial; ss=ss_initial; mu=mu_ending; iota=iota_ending; gamma_a=gamma_a_initial; gamma_b=gamma_b_initial; rho_a=rho_a_initial; rho_b=rho_b_initial; chi_b=chi_b_initial;
    else
        phi=phi_initial; theta=theta_initial; ss=ss_initial; mu=mu_initial; iota=iota_initial; gamma_a=gamma_a_ending; gamma_b=gamma_b_ending; rho_a=rho_a_ending; rho_b=rho_b_ending; chi_b=chi_b_ending;        
    end

    %--------------------------------------------------------------------------
    %  A) SOLVE THE GE EQUILIBRIUM FOR THOSE PARAMETERS
    %--------------------------------------------------------------------------

      guess=[0.9,0.1,0.02,0.001,11,2,1,6,1,2,0.01];
    geq_vec=@(input)geq_solver_ext(input,iota,mu,gamma_a,gamma_b,chi_b,rho_a,rho_b,phi,theta);
    [answer,value,exitflag] = fsolve(geq_vec, guess, options);
    if exitflag ~= 1
        disp('First-Order Optimality is not Achieved')
    end
    
    iah=answer(1);
    ial=answer(2);
    ibh=answer(3);
    ibl=answer(4);
    omegah=answer(5);
    omegal=answer(6);
    a=answer(7);
    v1h=answer(8);
    v1l=answer(9);
    nsa=answer(10);
    nsb=answer(11);


    %--------------------------------------------------------------------------
    %  B) RECOMPUTE ALL “INTERMEDIATE” QUANTITIES EXACTLY AS IN MAIN PROGRAM
    %--------------------------------------------------------------------------

    omega_avg=alphah*omegah+alphal*omegal;
    nbah=alphah*(1-iah*Xa(omegah));
    nbal=alphal*(1-ial*Xa(omegal));
    nba=nbah+nbal;
    nbbh=alphah*(1-ibh*Xb(omegah));
    nbbl=alphal*(1-ibl*Xb(omegal));
    nbb=nbbh+nbbl;
    
    sigmaah=nbah/nba;
    sigmaal=nbal/nba;
    sigmabh=nbbh/nbb;
    sigmabl=nbbl/nbb;
    
    mba=phi*(nsa/nba)^nu;
    msa=phi*(nba/nsa)^(1-nu);
    mbb=phi*(nsb/nbb)^nu;
    msb=phi*(nbb/nsb)^(1-nu);
    
    auxh=(alphah*mh*a+alphal*ml)*(alphah*Q(1,1)+alphal*Q(2,1))/((alphah*mh+alphal*ml)*(alphah*Q(1,1)*a+alphal*Q(2,1)));
    auxl=(alphah*mh*a+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2))/((alphah*mh+alphal*ml)*(alphal*Q(2,2)+alphah*Q(1,2)*a));
    gh=1+auxh*gamma_a*(iah*Xa(omegah)+(1-iah*Xa(omegah))*mba)+auxh*gamma_b*(ibh*Xb(omegah)+(1-ibh*Xb(omegah))*mbb);
    gl=1+auxl*gamma_a*(ial*Xa(omegal)+(1-ial*Xa(omegal))*mba)+auxl*gamma_b*(ibl*Xb(omegal)+(1-ibl*Xb(omegal))*mbb);
    
    g_avg=gh*(alphah*Q(1,1)+alphal*Q(2,1)/a)/(alphah*Q(1,1)+alphal*Q(2,1));
    rm=beta/(g_avg^(epsilon*zeta/(lambda+zeta)));
    rt=1/rm-1+delta;
    w=lambda*((eta/rt)^(eta/(zeta+lambda)))*((alphah*mh+alphal*ml)*g_avg*zbar)^(zeta/(lambda+zeta));
    
    oa=(msa*theta*(sigmaah*v1h+sigmaal*v1l)*gamma_a)/(1-(1-msa*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));
    ob=(msb*theta*(sigmabh*v1h+sigmabl*v1l)*gamma_b)/(1-(1-msb*theta)*rm*delta_p*g_avg^(zeta/(lambda+zeta)));
    
    zlbar=(alphah*mh+alphal*ml)*zbar/(alphah*mh*a+alphal*ml);
    zhbar=a*zlbar;    
    
    rd_exp_ah=(1-ss)*chi_a*(zbar^(zeta/(zeta+lambda)))*(iah^(1+rho_a))/(1+rho_a);
    rd_exp_al=(1-ss)*chi_a*(zbar^(zeta/(zeta+lambda)))*(ial^(1+rho_a))/(1+rho_a);
    
    rd_exp_bh=(1-ss)*chi_b*(zbar^(zeta/(zeta+lambda)))*(ibh^(1+rho_b))/(1+rho_b);
    rd_exp_bl=(1-ss)*chi_b*(zbar^(zeta/(zeta+lambda)))*(ibl^(1+rho_b))/(1+rho_b);
    
    rd_exp_h=rd_exp_ah+rd_exp_bh;
    rd_exp_l=rd_exp_al+rd_exp_bl;
    
    basic_h_share=rd_exp_bh/rd_exp_h;
    basic_l_share=rd_exp_bl/rd_exp_l;
    
    basic_share=(alphah*rd_exp_bh+alphal*rd_exp_bl)/(alphah*rd_exp_h+alphal*rd_exp_l);
    
    zhr=(alphah*Q(1,1)*zhbar+alphal*Q(2,1)*zlbar)/(alphah*Q(1,1)+alphal*Q(2,1));
    zlr=(alphah*Q(1,2)*zhbar+alphal*Q(2,2)*zlbar)/(alphah*Q(1,2)+alphal*Q(2,2));
    
    Yh=mh*gh*zhr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));
    Yl=ml*gl*zlr*((eta/rt)^(eta/zeta))*((lambda/w)^(lambda/zeta));
    
    rd_Y_h=rd_exp_h/Yh;
    rd_Y_l=rd_exp_l/Yl;
    
    patent_count=alphah*(iah+ibh)+alphal*(ial+ibl);
    trans_share=(alphah*iah*(1-Xa(omegah))+alphal*ial*(1-Xa(omegal)))*(msa*(1-((1-msa)*delta_p)^10))/((1-(1-msa)*delta_p)*patent_count)+(alphah*ibh*(1-Xb(omegah))+alphal*ibl*(1-Xb(omegal)))*(msb*(1-((1-msb)*delta_p)^10))/((1-(1-msb)*delta_p)*patent_count);
    
    
    if v1h*gamma_a >= oa && v1l*gamma_a>=oa
        disp('Both types keep applied patents');
    elseif v1h*gamma_a > oa && v1l*gamma_a<oa
        disp('High type keeps applied patents; low type sells applied patents');
    else
        disp('Both types sell applied patents');
    end
    
    if v1h*gamma_b >= ob && v1l*gamma_b>=ob
        disp('Both types keep basic patents');
    elseif v1h*gamma_b > ob && v1l*gamma_b<ob
        disp('High type keeps basic patents; low type sells basic patents');
    else
        disp('Both types sell basic patents');
    end

    % Package the five moments:
    M1 = struct();
    M1.basic_share=basic_share;
    M1.scope_avg = omega_avg;
    M1.RD_H    = rd_Y_h;
    M1.RD_L    = rd_Y_l;
    M1.trans   = trans_share;
    M1.growth  = g_avg^(zeta/(zeta+lambda))-1;

    % Compute percentage deviations (or percentage‐point for growth):
    deviations(i,1) = (M1.basic_share - M0.basic_share) / (basic_target2-basic_target) * 100;   % %Δ scope_H
    deviations(i,2) = (M1.scope_avg - M0.scope_avg) / (range_avg2-range_avg) * 100;   % %Δ scope_H
    deviations(i,3) = (M1.RD_H    - M0.RD_H   ) / (rd_sales_h2-rd_sales_h)    * 100;   % %Δ RD_H
    deviations(i,4) = (M1.RD_L    - M0.RD_L   ) / (rd_sales_l2-rd_sales_l)     * 100;   % %Δ RD_L
    deviations(i,5) = (M1.trans   - M0.trans  ) / (trans_target2-trans_target)   * 100;   % %Δ trans‐share
    deviations(i,6) = (M1.growth  - M0.growth ) * 100;                % Δ growth (pp)
end


%--------------------------------------------------------------------------
% 4) Generate Table 10
%--------------------------------------------------------------------------

rowNames = {
    'Patent Market (ϕ,θ)';
    'Matching Efficiency (ϕ)';
    'Sellers Bargaining Power (θ)';
    'R&D Tax Credit (σ)';
    'Production Cost (µ,ι)';
    'Rare Good Ideas (χ^a,χ^b,γ^a,γ^a,ρ^a,ρ^b)'
};

varNames = {
    'Basic Research'
    'Prod. Scope';
    'R&D(H)';
    'R&D(L)';
    'Patent Trade';
    'Growth'
};

Table10 = array2table(deviations, ...
    'RowNames',    rowNames, ...
    'VariableNames', varNames);

disp(Table10);
writetable(Table10,  "C:\Users\Yueyuan Ma\Box\Specialization\Replication_Files\Tables\Table10.xlsx", 'WriteRowNames', true);



